Fluid-Fluid and Fluid-Solid transitions in the Kern-Frenkel model from 
Barker-Henderson thermodynamic perturbation theory 
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^ ■ Abstract 

\ We study the Kern-Frenkel model for patchy colloids using Barker-Henderson second-order thermody- 

namic perturbation theory. The model describes a fluid where hard sphere particles are decorated with one 



■ patch, so that they interact via a square-well (SW) potential if they are sufficiently close one another, and 

O" 



if patches on each particle are properly aligned. Both the gas-liquid and fluid-solid phase coexistences are 
computed and contrasted against corresponding Monte-Carlo simulations results. We find that the pertur- 
bation theory describes rather accurately numerical simulations all the way from a fully covered square- well 
potential down to the Janus limit (half coverage). In the region where numerical data are not available 
(from Janus to hard-spheres), the method provides estimates of the location of the critical lines that could 
serve as a guideline for further efficient numerical work at these low coverages. A comparison with other 
techniques, such as integral equation theory, highlights the important aspect of this methodology in the 
present context. 
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I. INTRODUCTION 



Perturbation theory has a long and venerable history in the context of fluids and a detailed 
description of several different techniques is presented in classic textbooks 

HQ 

, and in excellent 

dedicated reviews 4]. 

Although the general idea dates back to a much earlier time, the first well established paradigm 
of first- and second-order perturbation theory was devised by Zwanzig Q] for simple fluids, later 
extended to polar fluids Qj . A similar analysis was carried out by Buff and Schindler in the context 
of solution theory 3]. 

In addition to these theories that assume the hard-spheres model as unperturbed system, other 
theories exist that rely on the van der Waals picture as a starting point, the best known of these 
being the Weeks-Chandler- Anderson (WCA) theory js-ll]. 



While the WCA theory has proven extremely powerful in many applications, for potential 
with hard-cores the original Zwanzig theory offers a natural scheme, hinging on an unambiguous 
potential separation. This was eventually put on firm ground by Barker and Henderson (BH) 



13], a rather unrealistic potential 



[3j, 0, [ijj who provided reliable estimates for square-well fluids 
in the framework of simple liquids, but much more physically sound when applied to the colloid 
domain. 

In the present paper, we will apply the BH thermodynamic perturbation theory to the Kern- 
Frenkel (KF) model for patchy colloids [l-jl . If]]. In this model 14], attractive circular patches are 
distributed on the surface of hard-spheres, and different spheres attract each other provided that 
any two patches on distinct spheres are suitably aligned, and the relative radial distance between 
meres is within the range of the attractive tail. 



the centers of the s 

While not new 16(, these systems have witnessed an impressive resurgence of interest in the 
last few years for a number of reasons. The first reason is due to the remarkable improvements in 
the chemical synthesis techniques that allows to decorate the surface of a colloid with great preci- 
sion and reliability, a feature that was not possible until few years ago. When combined with the 
additional advantage, as compared with their atomic counterpart, of an almost arbitrarily control 
of their size and interaction range, this makes patchy colloids very attractive for technolo gica l ap - 



plications, as elementary building blocks for self-assembly materials of the new generation 
An additional important reason hinges on the fact that patchy colloids may serve as a paradigm for 
systems with low valence, strong anisotropy, and highly directional interactions between particles, 



a feature that is common to many different systems, globular proteins being a notable example, 
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where the heterogeneity of the surface groups cannot be neglected even at the minimal level. 

Several examples of applications and improvements of the original BH scheme have been of- 
fered over the years. Verlet and Weiss discussed a comparison with numerical simulations and 
experimental results both for simple 



n 



191 ] and polar [20f| fluids; Gubbins and Gray 



perturbation scheme for molecular fluids; Chang and Sandler 



proposed a 



22j exploited it to develop an ana- 
lytical approximation for the square- well fluid valid within a particular interval of well amplitude; 
Zhang et al, 



24j used it to 



231 ] applied it to a square-well chain fluid, whereas Rotemberg et al. 
study the phase behavior of mixtures of colloidal particles and interacting polymers. More recently, 
Zhou [25( derived a simple procedure hinging on the BH scheme to locate the fluid-solid coexistence 



phase for a hard-core attractive Yukawa fluid, and Kalyuzhnyi et al. [26j, |27] tackled the single 
and multiple patchy colloids, similar to those treated in the present work, using a generalization 
of Wertheim's thermodynamic perturbation theory 



28H31[. 



The present work builds upon the methodology outlined in Rcf. 



321 ] to show that BH second- 



order perturbation theory can be successfully applied to patchy colloids, as represented by the 
Kern-Frenkel model Jj|. Besides thermodynamic quantities, such as virial equation of state and 
chemical potentials, the method allows a rather precise location of the fluid-fluid and the fluid-solid 
coexistence lines, in principle for arbitrary number and size of the patches. In this respect, the 
method competes in accuracy with integral equation theory on the same system without 
suffering from the unavoidable instabilities present in that case for low surface coverages and 
temperatures. This will be demonstrated by an explicit comparison with numerical simulations 
carried out [3314361] on the same system. 

The outline of the paper is as follows. In Section |H] the model is defined and in Section ITTT1 the 
used perturbation technique is described. Some technical details of the calculations are included in 
Appendices [A] and |Bj Section IIVI includes the method to compute the coexistence curves from the 
analytical results, with details of the numerical procedure included in AppendixO Section|V]briefly 
summarize some details of the Monte Carlo calculations, and Section IVTI includes descriptions of 
all results. Finally, Section IVIII summarizes the paper and provides some future perspectives. 



II. THE KERN-FRENKEL MODEL 



Consider a fluid formed by N particles in a volume V at temperature T, and assume that they 
can be described by the Kern-Frenkel model [jjj in its one-patch version (see Fig. [T]), where the 
orientation of the patch on each surface sphere 1 and 2 is identified by unit vectors rii and n2, 
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whereas the direction connecting centers of spheres 1 and 2 are identified by unit vector ri 2 . 

Two spheres of diameter a attract via a square-well potential of width (A — l)a and depth e if 
the directions of the patches on each sphere are within a solid angle defined by 6q, and repel each 
other as hard spheres otherwise. The pair potential has the form 

$(12) = (n 2 ) + $i (12) > ( 2J ) 
where the first term is the hard-sphere contribution 



^ (r)={ ' , (2.2) 




and the second term 

$i(12) = </»sw(ri 2 )*(ni,n 2 ,f 12 ) (2.3) 

is the orientation-dependent attractive part which can be factorized into an isotropic square-well 
tail 

{— e, a < r < Xa 
(2-4) 
0, Act < r 

multiplied by an angular dependent factor 

1, if ni • ?i2 > cos Qq and — n 2 • ?i2 > cos 6 
*(ni,n 2 ,r 12 ) = ^ . (2.5) 

I 0, otherwise 

The unit vectors nj(wj), (i = 1, 2), are defined by the spherical angles = (#j, ipi) in an arbitrarily 
oriented coordinate frame and ri 2 (Q) is identified by the spherical angle £1 in the same frame. 
Reduced units, for temperature T* = ksT/e, pressure P* = fiP/p and density p* = pa 3 , will be 
used throughout, with ks being the Boltzmann constant. For future reference, we also introduce 
the packing fraction r/ = irp* /6. Two particles then attract if they are within the range of the 
square-well potential and if their attractive surfaces are properly aligned with each other, and 
repel as hard spheres otherwise. 

The relative ratio between attractive and total surfaces is the coverage x that is related to the 
semi-angular width #o of the patch. This can be obtained as 

X 2 = (* (Ai, n 2 , r 12 )>^ 2 = (* 2 (n x , n 2 , r^))^ = sin 4 (|) , (2.6) 

where we have introduced the angular average 

(...). = tJ^-- (2 - 7) 
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III. BERKER-HENDERSON PERTURBATION THEORY 



The Kern-Frenkel potential defined in Eqs. (|2TT| . (12"^) . (1231) . (l?^) . and ([23]), leads to the 
natural separation into a reference one (the hard-sphere contribution) and an interaction term (the 
remaining, angular dependent, part) that is usually requested by the standard perturbation theory 
prescription [l[ [J . 

The presence of the hard-sphere potential for the reference part further suggests the Barker- 
Henderson (BH) scheme [l2] as the most suitable one for the present model. This has also the 
additional advantage that the free energy Fq for the reference system can then be computed in 
several ways, as further discussed below. 



The original scheme, due to Zwanzig 



5J, provided the first and second-order terms within the 



canonical ensemble, in the form of a high temperature expansion 

P(F-F ) m 3F 2 

where the first term is proportional to 1/T* the second to (1/T*) 2 . 

nnn 

Although formally correct, it was noticed by Barker and Henderson |3j, |4|, LL2| that the corre- 
sponding expressions were not useful for finite systems and a grand canonical ensemble derivation 
provided a much more convenient framework, where the results for the canonical ensemble could 
be eventually obtained by a Legendre transformation. 

To the best of our knowldege, the details of the computation for the second-order term were 
presented in Ref. only for isotropic potentials. As its generalization to angular dependent 
potentials proves to be instructive, we have outlined in Appendix lAl 

The first term poses no problem and is computed in Eq. (|A18p . When the perturbation param- 
eter 7 = 1 and particularized to the Kern-Frenkel potential given in Eqs. (|2.ip . (12. 2p . ()2.3fl . (|27 
and (I2.5P it becomes 

0Fi 1277 [ X ° 



drr'g Q {r)4>w{r)m {12))^ 2 . (3.2) 

Note that this term is negative because so is 4>sw(f)- 

The second term is much more involved, but one can apply the same procedure as the isotropic 



case 



32l | . as detailed in Appendix [Aj The result for the second term is reported in Eq. (]A19p . As 



in the isotropic case, however, this derivation is of little practical use in view of the presence of the 
three and four point distribution functions [3]. Barker and Handerson jl^ . devised then a simpler 
procedure to compute this term, based on a discrete representation of the radial axis distributions. 
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Again, the original procedure for spherically symmetric potentials is extended to angular dependent 
potentials in Appendix IB] 

The result for the second-order term is reported in Eq, (|B14p . In case of the Kern-Frenkel 
potential, it yields 



(3F2 6r] ( drj 



Xrr 



drr 2 g (r)4 w (r)([^(12)} 2 ) , (3.3) 



N a 3 \dP* 

where Pq = 0Pn/p is the reduced pressure of the HS reference system in the Carnahan-Starling 



approximation 



37]. 



This result is identical to that reported in Ref. [32 1 for a different radial part and it is known 

□ 

as macroscopic compressibility approximation [12|. Although the results given in Eqs. (|3.2p (first 
order) and (13. 3D (second order) are somewhat intuitive, being the natural extensions of the isotropic 
counterpart, a detail analysis of their derivations is important as it might help to improve a 
drawback of the method that will be discussed at the end of Section IIV1 



IV. FLUID-FLUID AND FLUID-SOLID COEXISTENCE CURVES 



Once the reduced free energy per particle /3F/N is known, all thermodynamic properties can be 
derived. In particular, the pressures and the chemical potentials can be derived from the standard 



thermodynamic identities 



PP 
P 

Pp 



d_(PF_\ 
V dri\N J 
d_ ( PF\ 
drf y~NJ 



(4.1) 
(4.2) 



The gas-liquid (fluid-fluid) and fluid-solid coexistence curves are determined by the equality 
of the temperature, pressure and chemical potential in the two coexisting phases. Since the two 
phases are in contact, the condition on the equality of the temperature is always fulfilled. Thus, at 
fixed temperature T* , we are left with the two conditions on the pressure and chemical potential. 

For the gas-liquid coexistence, the conditions are 



p;{T*, P * g ) = ir(T*,p*) 

M*(r\ P *) = M?(r*,pf) , 



(4.3) 
(4.4) 



where subscripts g, I indicate that the quantity is computed in the gas and liquid phase respectively. 
The solution of this system of non- linear equation gives p* = p*(T*) for the gas coexistence branch, 
and p\ = pf(T*) for the liquid coexistence branch. The hard-sphere reference part of the free energy 



(in excess with respect to the ideal gas) is assumed to be described by the Carnahan- Starling 
relation 

4r/ 2 - 3rj 3 

liquid 



N 



77(1-77)' 



(4.5) 



For the hard-sphere radial distribution function go (r) part 



the Verlet-Weis 



19 



39 



ap 



searing in Eqs. (|3.2p and (|3.3|) 



40l | is exploited. The details of the 



corrected Percus-Yevick solution 
numerical procedure are reported in Appendix [Cl 

A similar method can be applied to the fluid-solid transition, resulting into the conditions 



, P * f ) = p;(t*, p * s ) 

,p* f ) = p* s (T*,p* s ) . 



(4.6) 
(4.7) 



All calculations assume that the solid phase retains the crystal structure of the reference system, 
namely the face-centered cubic (fee) lattice. It is possible, especially at low T or low x where 
anisotropy effects are more relevant, that fee is not the most stable solid for the model; our 
coexistence results are still valid, although possibly relating to a metastable solid phase. We used 
Wood's equation 



41[ 



N 



2.1306 + 3 In 



V 



1 



+ ln 



pA 3 



(4- 



solid \--v/Vc P J V V 

for the solid free energy of the reference hard-sphere part, where r] cp = 7r\/2/6 is the fee volume 
fraction for closed packing. For go (r) in the solid phase, we use the orientation-averaged pair 



distribution function of Kincaid and Weis 



42] 



As a double check of the reliability of the numerical solution of Eqs. (|4.6p and (|4.7p . the critical 
points were also computed using the alternative, and more direct method, as a maxima of the 
Helmholtz free energy, that is from the system 



d 2 ($F_ 
dp* 2 \~W 
d 3 ( (3F 

w 



dp 




0. 



(4.9) 



(4.10) 



The solution provides T* and p* and are consistent with previous results, as they lay exactly on 
the top of the coexistence curves. 



8 



V. MONTE CARLO SIMULATIONS 



Standard Monte Carlo (MC) simulations in the NPT and in the grand-canonical (GC) ensem- 
bles have been implemented to evaluate the equation of state and the density dependence of the 
chemical potential for both the Janus and the SW model. Translational and rotational moves 
consist of random translation of ±0.1<7 and random rotation of ±0.1 rad of a randomly selected 
particle. In the GC study, insertion and deletion moves have been attempted, in average, every 
500 displacement moves. In NPT simulations, N = 500 particles were investigated. In GC sim- 
ulations box sizes were selected in such a way that the number of particles in the simulation box 
was would always larger than 500. Fluid-fcc coexistence lines were calculated via Gibbs-Duhem 
integration 



methods 



431 ] . starting from initial coexistence values at T* = 2 established via direct coexistence 
44| . Since at infinite temperature the KF model reduces to the hard sphere model, coex- 



istence pressures at T* = 2, a very high value for the K 
of the known HS values. We refer the reader to Refs. 44J, 



model, were searched for in the vicinity 



451 ] for the details of the procedures. We 



point out that all NPT simulations of the fee solid were carried out in a cubic box to constrain 
the system to retain the fee arrangement also in cases where the preferred structure would be a 
different one, possibly other lattices or a distorted fee. This choice was made to properly compare 
simulation results with the perturbation theory that assumes the cubic fee arrangement of the 
reference SW system. 



VI. RESULTS 



A. Equation of state and chemical potential 



In order to assess the performance of perturbation theory, we first compare results for pressure 
and chemical potential as derived from the BH scheme outlined in Sect Mil with numerical simula- 

from integral 



tions 
equation 
by Lado 



3314361] . These values were further compared with those derived in Ref . 



33 



IE ) th eory within the reference hypernetted chain (RHNC), following the method devised 



4614481]. In the square- well case, integral equation values were taken from Ref. [491] . 



The results are shown in Figures [2] (pressure) and [3] (chemical potential), for two representative 
values of the coverages, namely the square-well (x = 1) (top panels) and the Janus (x = 0.5) 
(bottom panels). In all cases, a value of A = 1.5 for the total extension of the well (in units of the 
hard-spheres diameter), was selected in order to compare with previous results. 

In the square- well case (top panels), the first selected temperature ksT/e = 1.4 corresponds 
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to a tem per ature above the critical temperature, while the last one ksT/e = 1.0 is well below 
(see Ref. [491] and references therein). In both cases, the performance of the BH thermodynamic 
perturbation theory is remarkably good, in agreement with previous findings on the square-well 
potential The dip in the curve for ksT/e = 1.0 indeed corresponds to the van der Waals 

loop typical of the coexistence region. In the case of chemical potential (FigE]), the ideal gas low 
density solution /3/i = ln(pcr 3 ) is also reported for comparison. Interestingly, while both Monte 
Carlo (MC) simulations and RHNC integral equation theory (IE) converge to the correct limit, 
the BH perturbation theory appears to underestimate the chemical potential in the whole range 
of densities. On the other hand, it provides the same quality results for all temperatures, even in 
those regions where integral equation theory are known to experience difficulties. 

Slightly less satisfactory results are obtained in the case of a Janus fluid, as shown in the bottom 



panels of both Figs. [2]and[3l Here the two 
are both above the critical temperature 



imiting temperatures fc^T/e = 0.9 and ksT/e = 0.55 
361 ] , as apparent from the absence of any loop in both 
the pressure and the chemical potential. The Janus phase diagram, however, is known to be 
anomalous [35J] , as a result of a competition with a micelle formation process that destabilizes the 



condensation one [36[]. In this case the BH thermodynamic perturbation theory (BH) does not show 
a well defined pattern as it overestimates the pressure for both temperatures (Fig[2] bottom panel), 
as well as the chemical potential for ksT/e = 0.55, but underestimates it for the higher temperature 
ksT/e = 0.9 (Fig. [3]bottom panel). While it is known that the BH compressibility approximation 
can be expected to display different performance at different densities due the presence of higher- 
order terms [3], the above behavior is more likely to be attributed to the anomalous behavior of the 
Janus phase diagram that perturbation theory cannot capture at the present stage. In spite of this, 
the performance of BH thermodynamic perturbation theory remains remarkable, especially in view 
of the difficulties experienced by integral equation theories at such low temperatures associated 
with low surface coverages. 

B. The fluid-fluid coexistence 

A very stringent test on the reliability of BH thermodynamic perturbation theory stems from 
the calculation of the fluid-fluid (gas-liquid) coexistence curves. This is depicted in Figure H] where 
the coexistence curves are computed by BH thermodynamic perturbation theory (solid lines), and 



contrasted with results from Monte Carlo numerical simulations (points), from Ref. [3g]. The 



considered coverages range from \ = 1.0, corresponding to the SW fluid, to \ = 0.5, corresponding 
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361 ] . As before, A = 1.5 was assumed 



to the Janus limit, and are identical to those considered in Ref. 
in all calculations. In the BH thermodynamic perturbation theory, further coverages down to the 
hard-sphere limit were also computed. In all cases, the critical points stemming from the BH 
calculations are also displayed as solid circles on the binodals. 

The performance of the BH method appears to be remarkably good. Both the vapor and the 
liquid branches of the numerical simulations are closely followed by the BH calculations, with an 



y 



exception of the Janus 



351 ] . as remarked. This is 



accuracy almost independent of the considered coverage, with the on 
case (x = 0.5) that is however known to have an anomalous behavior 
only apparently in contrast with results from chemical potential, reported in previous Section TVI Al 
where the BH results for chemical potential in the Janus case appeared to be less precise than in the 
SW case. On the one hand, a closer inspection reveals that BH results for each single coverage do 
indeed show a small quantitative discrepancy with the corresponding MC simulation, more or less 
uniform in the entire density-temperature plane. On the other hand, this latter feature constitutes 
an advantage in the method as a numerical solution of Eq. (|4.4p may provide accurate coexistence 
lines if both the vapor and the liquid chemical potentials have similar inaccuracies. This results 
is, nonetheless, comparable in accuracy with those stemming from reference hypernetted chain 



(RHNC) integral equation theory [33|, [34], |49( , with the additional advantage of a less computational 



and algorithmical complexity and, more importantly, of being able to access the critical region, 
including the critical point, that is one of the main shortcomings common to all integral equation 
theories. 

It is worth noticing how BH perturbation theory can provide an accurate prediction of the 
location of the coexistence lines even below the Janus limit, that is for x < 0.5, where extensive 
numerical simulations are so-far suggesting the fluid-fluid transition to be inhibited by a micelliza- 
tion process [36|. This could be useful for a future more focussed numerical calculation within a 



limited region of the phase diagram where BH theory predicts coexistence to occur. 



C. The fluid-solid coexistence 



Let us now turn to the fluid-solid coexistence, a calculation that has not been carried out so 
far for this model by any method. As illustrated in Sec JIVI and below, BH perturbation scheme 
allows this analysis with an effort, both computational and algorithmical, comparable with that of 
the fluid-fluid case. 

In the isotropic SW case with A = 1.5, the reference point for this calculation are those obtained 
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as early as in 1980 by Young and Adler [50J. Using molecular dynamics (MD), they reported 
a detailed study of the different crystal structures (fee, hep, and bec) with the corresponding 
Helmholtz free energies, thus arguing that fee and hep were the most stable structures within the 



entire temperature-density plane. Additiona 
and Gast [51], Serrano-Ulan and Navascues 



numerical simulations were later performed by Marr 



531 ] essentially confirming this 



52j, and Kiselev et al. 

scenario. A very detailed study of the entire phase diagram of the SW, was carried out by Liu et 
al. [54]. 



In Figj5j we report results from BH thermodynamic perturbation theory (solid line) along with 
results from Young and Adler (circles). 

While at high temperatures all calculations agree, discrepancies start to appear on cooling away 
from the hard spheres limit. In particular, the plateau appearing in the solid branch of MD calcu- 
lations indicates a fcc-fcc (or fec-hep) transition that is not accounted for in BH calculations, that 
assumed fee structures all the way, although in principle it could be done. In the BH calculations, 
in particular, the difficulty arises from the stability of the numerical scheme used for the solution 
of Eqs.dMI) and (J377j) . 

For lower coverages, no previous calculations on the Kern-Frenkel model exist to compare with. 
Figj6] illustrate the coverage dependence of the fluid-solid coexistence lines as computed from MC 
simulations (points) and from BH thermodynamic perturbation theory (lines). As in the fluid-fluid 
case, MC simulations have been obtained up to the Janus fluid (x = 0.5), whereas BH theory 
provides results even below that limit. Simulations below the Janus limit could be done, but are 
computationally more demanding. 

As in the SW case, even for lower coverages one might expect a structural transition at a certain 
density. Even assuming fee to be the most stable structure, the range of the potential associated 
with the value A = 1.5 used here, allows a fcc-fcc transition between one fee with only nearest- 
neighbors bonded, and a more denser one with even the next-to-nearest-neighbors are bonded. 
This is associated with the jump in density that is present in some of the plots of FigJBJ 



VII. CONCLUSIONS 



In this paper we presented the first Barker-Henderson (BH) perturbative calculation for the 
one-patch Kern-Frenkel model, and compared with specialized MC simulations. The BH method 
hinges on a second-order thermodynamic perturbation theory in the inverse temperature, allowing 
the calculation of the Helmholtz free energy within this approximation, and hence, of all ther- 
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modynamic quantities, such as the density and temperature dependence of pressure and chemical 
potential. A numerical solution has then been implemented to infer the fluid-fluid coexistence line 
(binodal) from the equality of pressure and chemical potential in the vapor and liquid phase at a 
given temperature. A similar procedure also provides the fluid-solid transition. 

When compared with numerical simulations, the BH predictions are found to be extremely 
reliable in the entire phase diagram, and for all coverages from the isotropic SW potential to the 
lowest considered coverage {x = 0.1) very close to the hard-spheres limit. This constitutes one of the 
main advantages with respect to, in principle, superior and more accurate theoretical methodologies 
hinging on integral equation solutions, that are typically affected by the impossibility of accessing 
the critical region, and by the numerical instabilities occurring at low temperatures associated 
with low coverages. Even at the quantitative level, BH results were found to be competitive with 
integral equation theories, in agreement with previous results on the isotropic SW fluid. 

The performance of BH is particularly noteworthy for coverages below the Janus limit, that 
is for x < 0.5, the most challenging region for numerical simulations in view of the tendency for 
the particles to form single and multi-layer clusters always exposing the hard-sphere surface in the 
outer region in order to maximize favorable contacts. This mechanism competes and destabilizes 
the condensation process and the interpretation of numerical simulation results become more and 
more obscure in that region. As a result, a clear scenario suggested by numerical simulations in 
this region is still missing. A better understanding could in principle be favored by our BH results 
that provide a well defined and restricted region of the temperature-density plane where indication 
of possible coexistences could be sought for. 

While in the present paper the BH method has been applied to a single patch Kern-Prenkel 
potential, the method could potentially be extended to higher number of patches with no difficulties. 
As a matter of fact, this has already been done in Ref. 32|] for two-patch colloids with Yukawa 



interactions for the attractive part. An inspection of the relevant equations (|3.2[) and ()3.3j) . however, 
suggest the result to be identical to the one-patch case at the same coverage. This means that the 
BH method, in the present form, is not capable of distinguishing between one and two-patches, at 
the same coverage, a feature that , conversely, is accounted for in both numerical simulations 



33 



35, 



36] 



and integral equation theory 

behavior present in the Janus limit of the sin gle patch 
coverage (x = 0.5) of the double patches model 



34f |. In particular, it cannot then account for the anomalous 



351 ] and not present in the corresponding 



36J] . This is rather surprising in view of the fact that 



a similar method, based on a low-density virial expansion, applied to a companion problem, was able 
to distinguish between single and double patches, albeit with a rather poor estimate for the fluid- 
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fluid transition 55)]. A promising approach in this respect appears to be the perturbative scheme 



211 ] . who considered an expansion 



devised for molecular fluids by Gubbins, Gray and others [2, 
in powers of the anisotropic part of the potential, in a way akin to that discussed in Appendix 
lAl often supplemented by_ a Pade approximant to improve the convergence of the expansion, as 
proposed by Stell et al 



We plan to investigate this and other points in details in future work. 
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Appendix A: The second order perturbation theory 



As explained in Ref.[3j, the most correct way of developing a perturbation expansion is in the 
grand-canonical ensemble. Assume a general potential of the form 



U y (l,...,N) = U (l,...,N)+'yU I (l,...,N) 

= Yl $7 (0' ) = Y (ij) + 7 Y $I fa') ' 

i<j i<j i<j 



(Al) 



where Uq(1, . . . , N) = ^o(u) is the unperturbed part and Ui(l, . . . , N) = J2ij ^i(u) is the 
perturbation part. Here < 7 < 1 is used as perturbativ parameter, and each coordinate % includes 
both the coordinate rj and patch orientation hi, so that i = (rj,nj). Also, j3 = \/{ksT) denotes 
the inverse of the thermal energy. 

Introducing the following short-hand notation 

N 



Jl,...,N J 



n**<(---)>c 



i=l 



(A2) 



for the integration over all particle coordinates, the grand-canonical partition function 



E 



N\A™ Jim 



-0U 7 



(A3) 



N=0 " • T 

(here At is the de Broglie thermal wavelength, and f2 7 is the grand-potential) can then be used to 
obtain an expansion of the Helmholtz free energy 



'8F y 



7=0 



1 2 

H 7 

2! 1 



d 2 F 1 

dj 2 



(A4) 



7=0 
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as follows [jj. 

Taking the derivative of In Q 7 at fixed chemical potential fi, one has, using Eq. (|Aip 



[£hf ' 



= \ I 7rH*M 12 )]Mi2) , 



where 



-, +00 

Pl {i... h ) = —y: 



(iv-zoiA^A,..,* 



-/3C/ 7 



(A5) 



(A6) 



The second derivative is somewhat more laborious [3J, and one obtains after some algebra 



1 



9 



In Q 7 j = ^J i2 ^2 [-/^7 (12)] P 7 (12) + \J i2 [f^ H3* 7 ( 12 )] J /°7 (12) (A7) 
+ / |-[-/3$ 7 (12)]|-[-/3$ 7 (23)]p 7 (123) 

Jl,2,3 <?7 07 

+ 7 / #" i-^i (I 2 )] #" ["^7 (34)] [p 7 (1234) - p 7 (12) p 7 (34)] . 

4 .71,2,3,4 °1 °1 

The free energy F 7 is then obtained by considering 7 as an additional thermodynamical variable, 



and by performing the appropriate manipulations 



One then has 



F 1 = imN - k B TlnQ. 



7 ' 



(A8) 



and 



N = k B T 



d_ 



In Q 7 



(A9) 



where, for notational simplicity, here we do not distinguish between the canonical and grand- 
canonical number of particles N. Then 



k R T 



d_ 

d-f 



— In Q 7 



<9f2 7 \ ( dp\ ( dji 



9v y 7 \&yj \d P/ 



and hence, using the chain rule 



d l)p \ d v) p \dp) -y 



one gets 



k B T 



d_ 

O7 



— lnQ 7 



k B T 



d_ 

O7 



— In Q 7 



k R T 



d_ 



— In Q 7 



that, with the help of Eq. (IA9j) . leads to 



8E, 



-k B T 



d_ 



lnQ 7 



(A10) 



(AH) 



(A12) 



(A13) 
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where the right-hand-side is given by Eq. (|A5|) . 

For the second derivative, one proceeds as before, to obtain 



d 2 F 1 

dj 2 



d 2 



+ k B T 



( d 2 



Using Eq, (lA9j) and the relation 



dp \ d 



dps 



In Qr 



d_ 

dp 13 y dP J dp ' 



one finds 



2 



in a 



djdp 

Substituting in Eq. (|A14j ). one finds 



dp \ d 



8 2 F 1 
dj 2 



1 dP) dp 

N ( dp 

+ 



~[-^ 7 (12)]^(12) 



V 2 \dP 



d_ 

dp 



iA[_^ 7 (i2)]p 7 (12) 



(A14) 



(A15) 



(A16) 



, (A17) 



where the first term on the right-hand-side is given by Eq. (|A7p . 

The first and second order solutions, can be finally particularized to the potential form given 
in Eq.flXT]), so that Eqs.flAEJ) and (lAHfl) lead to 

d 



dj 



(/3F 7 ) 



7=0 



l -pN j dr 12 (/3$i(12)) wl|Wa g (12) 



(A18) 



and Eqs.flSH) and flSE]) leads to 

d 2 



d 7 2 (W 



7=0 



l -Np [ dr 12 ([-pa* (12)] 2 \ go (12) 

Np 2 j dr 12 dr 13 (12)] [-(3^ (23)])^^ g (123) 

- A Np* j dr 12 dr 13 dr M (12)] [-/3$ t (34)])^ 

^p 2 J dr 12 ($i(12)) w g (12) 



2 



(A19) 



[go (1234) -go (12) g (34)] 



Appendix B: The Barker-Henderson discrete representation 



As in the spherically potential case, the above expressions are, however, not very useful for 
practical computation, due to the high complexity involved in the calculations of the three <7o(123) 
and four go(1234) point correlation functions. 

Following the original work by Barker and Henderson, we return to the canonical partition 
function 

1 



Q 



N\A 3 T N /,. 



-PU(X,...,N) 



.N 



N\A 3 T N 



Z 



-PF 



(Bl) 
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that is related to the configurational partition function Z and to the Helmholtz free energy F. 
The intermolecular distance axis rij is divided in intervals (0, r\), (r±, r 2 ), ■ ■ ■ , (77, 77+1), ... in such 
a way that there are Ni distances in the I— th interval (77, 77+1 )• The total potential U appearing in 
Eq. (jBl|) can then be written as a sum over the different intervals with the respective multiplicity 

U(1,...,N) = N{§ (r h {Q, u;},) , (B2) 
i 

where <&(ri,{Q,u}}i) is the average potential in the I— th interval (assumed to be constant), and 
uj}[ are the set of orientational angles included in the same interval. 
Again we assume that the potential can be split into a hard-sphere part plus a tail 

¥ (n, {n, uj} t ) = 0o fa) + ¥1 ( n , {n, w },) . (B3) 

Introducing the average over the unperturbed system having Zq as configurational partition 
function 

(•••)o = 4- E f dr 1 ---dr N e-PZ l KiUn) ) (B4 ) 
Z ° n u n 2 ,... Jr 

where the symbol R indicates that the integral is restricted to configurations having JVj intermolec- 
ular distances in the interval (77,77+1), the Helmholtz free energy can be written in terms of that 
of hard-spheres Fq as 

(3F = pFo-hx^e-^NiMnWrti)^ J . (B5) 

Note that the angular average over the {0} variables is included in the average (IB4D . 
Use of the cumulant expansion 

- In (e~ Xx ^ = A (x) — (<x 2 > - (x) 2 ) + . . . (B6) 

leads to 

p{F- Fq) = PF X + PF 2 + . . . , (B7) 

where 

p Fl = E((^i(n,{^,^)) M ) o , (B8) 

and where 

PF 2 = -5E((W5 1 (r, ) {n )ll ;} l )iSi I (r m ,{fi )ll ;}J} M ) o . (B9) 

lm 
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As [lj, 

rn+i 



(iV z ) = 2vrpiV / drr'g (r) , (BIO) 
the first order term becomes 

m = l - P N f dr^(r)03$i(r,n,a*,^)>^ (Bll) 
that, of course, coincides with Eq. ()A18j l. 



For the second term (|B9p . an approximation is required as the effect of three and four-body 
interactions is included. Following Ref . \}\ • we assume molecules in different intervals to be uncor- 
rected 

(JV,# m ) - (JV,) (N m ) = l^m, (B12) 

and the fluctuations within a given interval, being related to the hard-spheres compressibility 

(N?) -ml = (N^kaT^j (B13) 
Substitution of Eqs. ([EI]) and (lBl~3~]) into Eq.dHH]), along with Eq. flBTO]) . leads to 

m = -IkaTpN^j J dr 50 (r)([/3<I>i(r,O, Wl , W2 )] 2 )^^ , (B14) 



which is the extension of the Barker-Henderson result 



121 ] to angular dependent potentials. 



Appendix C: Determination of the phase coexistence curves 

To illustrate how the phase coexistence curves are found numerically, we consider in the following 
the phase separation into a gas and a liquid phase; the fluid-solid coexistence curve is determined 
correspondingly. Our aim is to solve Eqs. (|4.3j) and (|4.4p for the two unknown particle densities p* 
and Pi of the gaseous and liquidus phase, respectively. Using the common tangent construction, 
the concentration of the density of the gaseous and liquidus phase can be found geometrically 



561 ] . In practice, however, p* and p\ is determined numerically by solving Eqs. (|4.3p and (|4.4p 
simultaneously using a nonlinear root finding algorithm. To illustrate this procedure, we rewrite 
Eqs. (03]) and (JOJ as 

fci (p* g ,pt) = P* g {T*,p* g ) - If (T*,pl) = (Cl) 
h2 {p* g , Pi) = Pi (T*, p* g ) - i4 (T*, p\) = , (C2) 
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where we have introduced the functions h\ (p*,p*) and \i2 \p* g ,p\). Since T* is kept fixed in the 
following, we have written h\ and /12 as function of p* and p\ only. By defining x = (p*,p^) t 
and h = (hi^hz)*, where the subscript t denotes the transposed matrix, our task of finding the 
concentrations of the two coexisting phases at constant T* is expressed in the form, 

h(x) = 0. (C3) 

This set of two nonlinear integral equation with two unknown variables is solved by using a well- 



tested implementation of the Newton-Rap hson method |57[, which solves Eq. (|Q3[) iteratively as 
briefly described in the following. First, a start value xq is chosen, and the gradient V/i(xo) is 
calculated. The new value x\ is found by a downhill step, 

x\ = xq — J~ 1 /i(x ) . (C4) 

Here, J is the Jacobian matrix which incorporates the partial derivatives of h\ and hi- This step 
is repeated, x\ — > x*2 — > £3 — > • • ., until the fix point x n = x* with 

h(x*) = 6, (C5) 

is found. It is important to note here that the root finding procedure requires the evaluation of h(x) 
at discrete points xi only. The nonlinear solver just steps down h{x) until Eq. (IC3P is fulfilled to 
a prescribed threshold. Since the evaluation of h(x) at x = X{ demands the calculation of several 
integrals, see Eqs. (|3.2p and ()3.3p . h(x) cannot be expressed in an analytical form. Hence, the 
nonlinear solver calls a subroutine which calculates both the free energy and its gradient for each 
iteration step X{. The free energy is evaluated using the Chebyshev quadrature and the derivatives 
in Eq. ()C4|) are calculated using Ridder's implementation of Neville's algorithm {5^ . 

After having found the two coexisting densities p* and p^ at a given T* , this procedure is 
repeated for a set of temperatures to map out the gas-liquid coexisting curve. The fluid-solid curve 
is calculated in exactly the same manner by equating the chemical potential and the pressure of 
the fluid and solid phase, Eqs. (|4.6|) and (|4.7|) . respectively. 
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FIG. 1: The Kern-Frenkel potential in the case of a single patch. The surface of each sphere is partionated 
into an attractive part (color code: green) and a repulsive part (color code: red). Units vectors rii and 112 
identify the directions of each patch, whereas the unit vector ?i2 join the centers of the two spheres, directed 
from sphere 1 to sphere 2. The particular case shown corresponds to a 50% fraction of attractive surface 
(coverage x = 0.5). 
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FIG. 2: Reduced pressure j3P/ p as a function of reduced density per 3 in the case of a square- well fluid with 
coverage x = 1 (top panel), and in the case of a Janus fluid with coverage \ = 0.5 (bottom panel). A value 
of A = 1.5 is used. Results from BH thermodynamic perturbation theory (BH) are compared with Monte 
Carlo simulation (MC) and with RHNC integral equation theory (IE). Different curves refer to different 
temperatures as shown. 
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FIG. 3: Reduced chemical potential as a function of reduced density per 3 in the case of a square-well 
fluid with coverage x = 1 (top panel), and in the case of a Janus fluid with coverage x = 0.5 (bottom panel). 
A value of A = 1.5 is used. Results from BH thermodynamic perturbation theory (BH) are compared 
with Monte Carlo simulation (MC) and with RHNC integral equation theory (IE). Different curves refer to 
different temperatures as indicated. The low-density ideal gas limit (light dashed line) is also depicted. 
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FIG. 4: The fluid-fluid coexistence curves as computed from BH perturbation theory and compared against 
numerical simulations. Lines are from perturbation theory, points from numerical simulations, for A = 1.5 
from Ref. [36(. All coverages from \ = 1-0 (SW case) to x = 0-0 (HS case) are depicted in the former case, 



whereas numerical simulations are in the range 0.5 < x < 1.0, that is from the Janus to the SW limit. 
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FIG. 5: Fluid-solid coexistence for the case of the SW potential (x — 1.0) with A = 1.5. Results from 
Barker-Henderson (BH) perturbation theory are contrasted with molecular dynamics (MD) data by Young 
and Adler 
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FIG. 6: Coverage dependence of the fluid-solid coexistence curves. Again A = f .5 and considered coverag 
are from x = 1-0 (SW case) to x = 0.1 for Barker- Henderson perturbation theory (lines) and from X = 
to x = 0.5 (Janus) for Monte Carlo simulations. 



27 



